Vitamin D and thyroid function: A mendelian randomization study

Background Numerous organs, including the thyroid gland, depend on vitamin D to function normally. Insufficient levels of serum 25-hydroxyvitamin D [25(OH)D] are seen as a potential factor contributing to the emergence of several thyroid disorders, however, the causal relationship remains unclear. Here we use a Mendelian randomization (MR) approach to investigate the causal effect of serum 25(OH)D concentration on the indicators of thyroid function. Methods We conducted a two-sample MR analysis utilizing summary data from the most extensive genome-wide association studies (GWAS) of serum 25(OH)D concentration (n = 443,734 and 417,580), thyroid-stimulating hormone (TSH, n = 271,040), free thyroxine (fT4, n = 119,120), free triiodothyronine (fT3, n = 59,061), total triiodothyronine (TT3, n = 15,829), as well as thyroid peroxidase antibody levels and positivity (TPOAb, n = 12,353 and n = 18,297), low TSH (n = 153,241), high TSH (n = 141,549), autoimmune hypothyroidism (n = 287,247) and autoimmune hyperthyroidism (n = 257,552). The primary analysis was conducted using the multiplicative random-effects inverse variance weighted (IVW) method. The weighted mode, weighted median, MR-Egger, MR-PRESSO, and Causal Analysis Using Summary Effect estimates (CAUSE) were used in the sensitivity analysis. Results The IVW, as well as MR Egger and CAUSE analysis, showed a suggestive causal effect of 25(OH)D concentration on high TSH. Each 1 SD increase in serum 25(OH)D concentration was associated with a 12% decrease in the risk of high TSH (p = 0.02). Additionally, in the MR Egger and CAUSE analysis, we found a suggestive causal effect of 25(OH)D concentration on autoimmune hypothyroidism. Specifically, each 1 SD increase in serum 25(OH)D concentration was associated with a 16.34% decrease in the risk of autoimmune hypothyroidism (p = 0.02). Conclusions Our results support a suggestive causal effect which was negative in direction across all methods used, meaning that higher genetically predicted vitamin D concentration possibly lowers the odds of having high TSH or autoimmune hypothyroidism. Other thyroid parameters were not causally influenced by vitamin D serum concentration.


Introduction
The human body receives vitamin D through synthesis in the skin after exposure to sunlight (a major source of vitamin D) or by the intake from food (only 5-10% of vitamin D is acquired in this way) [1].Vitamin D comes in five different forms (vitamin D1-D5).The two most important of these fat-soluble secosteroids for humans are vitamin D2 (ergocalciferol) and vitamin D3 (cholecalciferol).Calcitriol (also known as 1α,25-dihydroxyvitamin D (1,25(OH) 2D), 1α,25-dihydroxyvitamin D3 or 1,25-dihydroxycholecalciferol) is the active form of vitamin D3 produced in the liver from calcidiol (also known as 25-hydroxyvitamin D [25(OH)D] or calcifediol) by the action of the enzyme vitamin D 25-hydroxylase [2].Calcitriol has several important functions, such as immunosuppressive and anticancer effects [3,4], and it is involved in the regulation of phosphate and calcium levels, bone mineralization [5], neuromuscular and immune functions, regulation of cell growth, and the expression of more than 1000 genes [6,7].
The influence of diet on the functioning of the thyroid gland has been extensively documented in the literature.The well-known causal influence of iodine and selenium deficiency on thyroid dysfunction has motivated numerous studies on the effect of iodine and selenium supplementation on thyroid health [8][9][10][11], however, access to vitamin D (either from food or sunlight) is also crucial for the proper functioning of the thyroid gland.Meta-analyses [12,13] and observational studies (reviewed in Babić Leko et al. [14]) have shown that vitamin D insufficiency may be a risk factor for the development of hypothyroidism, autoimmune thyroid disorders and thyroid cancer.A large longitudinal case-control study based on electronic health records has found that vitamin D supplementation resulted in an overall decrease in TSH levels and lower rates of hypothyroidism detection [15].A recent randomized controlled trial (RCT) has found that vitamin D supplementation reduced the incidence of hypothyroidism in females, but not in males [16].Moreover, RCTs examining the impact of vitamin D supplementation on thyroid function in those suffering from autoimmune thyroid diseases found a decrease in anti-thyroid antibodies following vitamin D supplementation (reviewed in [14]).
The relationship between vitamin D and thyroid, however, is still poorly understood in many respects.Namely, due to the many inconsistencies between the studies (reviewed in [14]), it is still not clear how vitamin D affects thyroid hormones, thyroid stimulating hormone (TSH) and anti-thyroid antibody levels.Prior to implementing these findings in a clinical setting, it is essential to conduct additional research into the causal relationships that underpin these associations.The use of different assays to measure serum 25(OH)D levels and the possible confounding effects of age, sex, dietary practices, body mass index (BMI), smoking habits, and the season during which samples were taken are among the variables that could significantly contribute to discrepancies between the studies.
An approach that can further exclude the effect of the confounding factors, is Mendelian randomization (MR), which utilizes genetic variants as instrumental variables in order to estimate the causal effect of exposure on an outcome of interest.It is considered that this kind of methodology can mimic a randomized controlled trial [17,18].To date, only two studies using MR methodology investigated the association between vitamin D (mainly determined by serum calcidiol 25(OH)D levels) and thyroid function.MR methodology was used to investigate the association between serum 25(OH)D concentration and the risk of thyroid cancer, hypothyroidism, and hyperthyroidism [19], and serum 25(OH)D concentration and serum thyroid peroxidase antibodies (TPOAb) (the study included participants from China) [20].Additionally, a recent preprint [21] has investigated the association between serum 25(OH)D concentration and the risk of hypothyroidism, Hashimoto's thyroiditis and the levels of thyroid-stimulating hormone (TSH) and free thyroxine (fT4) using MR.However, to date, the association between serum 25(OH)D concentration and anti-thyroid antibody levels has not been investigated in populations of European ancestry.Furthermore, the associations between serum 25(OH)D concentration and the levels of free triiodothyronine (fT3) and total triiodothyronine (TT3) have not been investigated using MR methodology.Thus, the aim of this study was to assess the causal effect of serum 25(OH)D concentration on a comprehensive set of thyroid function parameters using MR methodology and the largest genome-wide association study (GWAS) summary statistics published to date.We investigated the following indicators of thyroid function; reference-range TSH, fT4, fT3, TT3, low TSH (levels below cohortspecific reference range), high TSH (levels above cohort-specific reference range), autoimmune hypothyroidism, autoimmune hyperthyroidism, TPOAb levels and TPOAb positivity.

Study design
We conducted a two-sample MR study using GWAS summary data to examine the possible causal effects of genetically predicted vitamin D levels on thyroid function and thyroid disorders.We leveraged extensive GWAS summary data for reference-range TSH, fT4, fT3, TT3, low TSH, high TSH, TPOAb, TPOAb positivity, autoimmune hypothyroidism and autoimmune hyperthyroidism (Table 1).Our methodological approach follows the one of our previous publication [22].For every causal relationship we examined, we used data from two independent but homogenous European ancestry populations.Our study encompassed a total of 20 MR analyses, aiming to shed light on the relationship between 25(OH)D concentration and thyroid function.Our research methods and reporting align with the STROBE-MR guidelines, as detailed in S3 Appendix.Our study did not need ethical approval given that our investigation used publicly available summary statistics.

Data Sources for Vitamin D
Two large GWAS on serum 25(OH)D concentration were published in 2020, a month apart from each other.Both of them utilized the UK Biobank data, leading to 417,580 participants in the Revez et al. [24] analysis and 443,734 participants in Manonusaki et al. [23], as the latter one also included an additional cohort of European descent participants.In order to exclude the bias relating to the choice of the GWAS, we have performed our analyses both with Revez et al. [24] and Manousaki et al. [23] GWAS.Both studies utilized a linear mixed model for the association analysis, however, the covariates and the transformations used differed significantly between the two studies (detailed in S1 Appendix).To estimate the level of agreement between the results of the two GWAS, we extracted the mutual instrumental SNPs and calculated the correlation coefficients for the beta estimates and p-values.Given that the summary association data for 25(OH)D concentration from both GWAS were missing rsIDs, we have imputed this information using the SumStatsRehab program [28].

Data sources for thyroid parameters
We utilized summary association data for reference-range TSH, fT4, fT3, TT3, low TSH and high TSH from the latest and most comprehensive GWAS meta-analysis on thyroid function from the ThyroidOmics Consortium [25].A comprehensive overview of the study cohorts and methodologies adopted can be found in Sterenborg et al. [25] and S1 Appendix.Summary association data for TPOAb levels and TPOAb positivity were obtained from an earlier Thyroi-dOmics consortium GWAS [26].
Autoimmune hypothyroidism and autoimmune hyperthyroidism GWAS summary statistics were obtained from the FinnGen consortium [27].The phenotypes used in this study were "hypothyroidism, strict autoimmune" and "autoimmune hyperthyroidism".FinnGen GWAS on hypothyroidism included 287,247 Finnish adults, of which 36,321 were cases and 250,926 were controls.FinnGen GWAS on hyperthyroidism included 257,552 Finnish adults, of which 1,621 were cases and 255,931 were controls.Participants classified as having autoimmune hypothyroidism included those who already received treatment, however, the same information was not present for participants in the autoimmune hyperthyroidism group.

Two-sample mendelian randomization analysis
In order to confirm the significance of the instrument's association with the exposure, we calculated the share of trait variance explained by the genetic instruments found (R 2 ) and evaluated the strength of these instruments using the F statistic.
The datasets for the exposure and outcome were aligned to ensure that the impact of identified single nucleotide polymorphisms (SNPs) on both the exposure and outcome was associated with the same effect allele.We used the Mendelian Randomization Pleiotropy RESidual Sum and Outlier (MR-PRESSO) technique for each combination of exposure and outcome to account for uncorrelated horizontal pleiotropy [29].After removing the SNPs that were shown to be horizontal pleiotropic outliers, we recalculated the causal relationship between the exposure and the outcome.
Our primary technique of analysis was a two-sample MR analysis using the multiplicative random effects inverse variance weighted (IVW) method [30].This method combines the influence of exposure on the outcome across various genetic variants and includes an overdispersion parameter into the variance of the (IVW) estimate, allowing the variance to expand in the presence of heterogeneity.While accommodating for heterogeneity among the SNP estimates, the multiplicative random effects IVW technique yields effect point estimates that are identical to those of the fixed effect IVW model [31].

Sensitivity analysis
For MR to produce valid causal inferences, the genetic instruments must satisfy three critical assumptions of instrumental variables: 1) Relevance, indicating a robust link between the instruments and the exposure; 2) Exchangeability, ensuring the instruments do not correlate with confounders that could affect both the exposure and the outcome; and 3) Exclusion, confirming that the instruments affect the outcome exclusively through their effect on the exposure [18].To verify our adherence to the relevance assumption, we computed the R 2 and the F-statistic for each exposure.Additionally, we used Cochran's Q statistic for each pair of exposure and outcome to check for heterogeneity in the SNP-specific estimates.
In addition to the primary multiplicative random-effects IVW approach, we validated our findings through MR analysis employing various alternative methods, each based on distinct model assumptions: 1) the weighted median 2) the weighted mode, and 3) the MR-Egger method.The weighted median method yields a reliable estimate when at least 50% of the weight is derived from valid instruments, making it resilient to outliers among the instruments [32].The weighted mode method operates on the assumption that the association estimate occurring most frequently is not influenced by pleiotropy, implying that it reflects the actual causal effect [33].We produced scatter plots to depict the impact of SNPs on the exposure in comparison to their impact on the outcome for each pair of exposure and outcome, utilizing the 'TwoSampleMR' package in R [34].Moreover, for each pair of exposure and outcome, we carried out the MR Steiger directionality test, which assesses whether the variance accounted for by the instrumental SNPs in the outcome is smaller than that in the exposure, to confirm the reliability of the directional analysis [35].We carried out this test using the 'directionali-ty_test' function within the 'TwoSampleMR' R package.
In addition to other methods, we utilized the Causal Analysis using Summary Effect Estimates (CAUSE) method [36], a Bayesian MR approach which, in contrast to traditional MR methods, utilizes the summary association estimates of all genome-wide variants available.CAUSE effectively addresses both uncorrelated horizontal pleiotropy (where a variant independently influences both outcome and exposure) and correlated horizontal pleiotropy (where a variant affects both the exposure and the outcome through a shared heritable factor).This approach makes the assumption that a certain percentage, denoted by q, of the variants exhibit correlated pleiotropy and that all other variants have uncorrelated pleiotropic effects.As a prior on q, this assumption is included in the model.Estimates of the posterior distribution for the sharing model and the causal model are generated using CAUSE.It then establishes whether the observed relationship between exposure and outcome is more likely the consequence of causality than correlated horizontal pleiotropy by determining whether the data fits the causal model more closely than the sharing model.Differing from other MR techniques, CAUSE utilizes complete summary association data from all variants, applying linkage disequilibrium pruning at at r 2 < 0.01 with P < 1 × 10⁻ 3 , thus enhancing the analytical power of MR analysis.
In all MR analyses, a two-sided p < 0.05 was considered statistically significant.We did not adjust for the family-wise error rate (FWER) in a traditional frequentist approach, as our selection was based on a specific subset of thyroid function parameters, and we further examined the presumed causal effects through a Bayesian framework.

Results
On average, the exposure analyses with the Manousaki et al.GWAS utilized 19.6 fewer instrumental SNPs (ranging from 27 to 89) than the exposure analyses with the Revez et al.GWAS (ranging from 38 to 113).The number of common instruments present in both exposure datasets was 33.Pearson's correlation coefficient between the beta estimates of these SNPs was equal to 0.44 (p = 0.01) while Kendall's correlation coefficient between the p-values of these SNPs was equal to 0.72 (5.58×10 −11 ).Because of the differences in methodological approaches and differences in obtained variant effects, we assumed that the subsequent MR analyses with these two exposure datasets would yield different results.

Serum 25(OH)D concentration from Manousaki et al. GWAS as exposure
We found a total of 109 independent SNPs for instrumental variables, which collectively accounted for 2.68% of the phenotypic variation (R 2 ) in 25(OH)D concentration, as detailed in S1 Table in S2 Appendix.Each chosen SNP had F-statistics above 10, with a median of 42.25 and a range of 34.60-64.27.Following the harmonization process, two SNPs were removed from all analyses due to being palindromic with intermediate allele frequencies, while an additional SNP was removed from the TPOAb analyses as noted in S2 The multiplicative random-effects IVW MR analysis indicated a causal link between genetically predicted serum 25(OH)D concentration and hyperthyroidism (S6 Table in S2 Appendix).Nonetheless, this causal relationship was not confirmed by sensitivity analyses or CAUSE analysis.The latter, by incorporating all genome-wide SNPs, showed that the causal model did not fit the data better than the sharing model (delta_elpd = 0.54 > 0, p = 0.78).

Serum 25(OH)D concentration from Revez et al. GWAS as exposure
We identified a total of 115 independent SNPs for instrumental variables, which collectively accounted for 3.52% of the phenotypic variation (R 2 ) in 25(OH)D concentration, as detailed in S1 Table in S2 Appendix.Each SNP chosen had F-statistics above 10, with a median of 43.13 and a range of 35.03-70.86.Following the harmonization process, two SNPs were removed from all analyses due to being palindromic with intermediate allele frequencies, while one additional SNP was removed from the high TSH analysis and two from TPOAb analyses as noted in S2 The multiplicative random-effects IVW MR analysis, as well as MR Egger analysis, showed that genetically predicted serum 25(OH)D concentration was causally associated with high TSH (Table 2 and Fig 1).Each 1 SD increase in serum 25(OH)D concentration was associated with a 12% decrease in the risk of high TSH (p = 0.0197) (Table 2 and S7  .The influence of genetically predicted 25(OH)D concentration on high TSH levels showed consistency in direction through both primary and sensitivity MR analyses, despite lacking statistical significance in the weighted median and mode analyses (S7 Table in S2 Appendix).Finally, the direction of the causal link was further validated by CAUSE analysis, which utilized all genome-wide SNPs to demonstrate that the causal model more accurately represented the data compared to the sharing model (delta_elpd = -0.86< 0).However, the disparity in fit did not achieve statistical significance (p = 0.31) (Table 2 and S1  Additionally, in the MR Egger analysis, genetically predicted 25(OH)D concentration was found to have a suggestive causal link with autoimmune hypothyroidism.Specifically, each 1 SD increase in serum 25(OH)D concentration was associated with a 16.34% decrease in the risk of autoimmune hypothyroidism (p = 0.02) (Table 2 and S7  However, there were indications of heterogeneity between variant estimates (Cochran's Q = 169.88,p = 3×10 −5 ) as well as uncorrelated horizontal pleiotropy (MR-PRESSO Global test p<0.0003),but MR-PRESSO was not able to identify any outliers (S5 Table in S2 Appendix).The effect of genetically predicted 25(OH)D concentration on autoimmune hypothyroidism showed uniformity in direction throughout both primary and sensitivity MR studies, yet failed to reach statistical significance in the multiplicative random-effects IVW, weighted mode, and weighted median analyses (S7 Table in S2 Appendix).
Finally, the direction of the causal link was additionally substantiated by CAUSE analysis.Utilizing all genome-wide SNPs, it indicated that the causal model aligned more closely with the data than the sharing model (delta_elpd = -0.56<0).Nonetheless, the difference in fit did not achieve statistical significance (p = 0.36) (Table 2 and

Discussion
Our results support a suggestive causal effect of genetically predicted serum 25(OH)D concentration on high TSH levels and autoimmune hypothyroidism.The results were noted in both frequentist and Bayesian analyses.Each 1 SD increase in genetically predicted serum 25(OH)D concentration was associated with as much as a 12% decrease in the risk of high TSH and a 16.34% decrease in the risk of autoimmune hypothyroidism.We found no causal effect of Consistent with existing literature, our findings contribute to the expanding research indicating a significant link between serum 25(OH)D concentration and thyroid function.Studies by Chailurkit et al. [37] and Kim et al. [38] found that higher 25(OH)D concentration was associated with lower TSH levels, particularly in younger individuals, and that low vitamin D status was associated with Hashimoto's thyroiditis.Similarly, Metwalley et al. [39] and Ma et al. [40] reported that vitamin D deficiency was prevalent in children and adolescents with autoimmune thyroiditis and that a lower 25(OH)D concentration was associated with an increased risk of developing autoimmune thyroid diseases.Additionally, our findings align with those of a large case-control study by Mirhosseini et al. [15], which reported that vitamin D supplementation led to fewer detected cases of hypothyroidism (including both clinical and subclinical forms) and an overall decrease in TSH levels.These findings suggest that maintaining adequate 25(OH)D concentration may play a role in reducing the risk of high TSH levels and autoimmune hypothyroidism.This growing body of observational evidence, paired with our causal results, could provide informed decisions in a clinical setting where vitamin D supplementation can be used as a cost-effective and safe strategy to improve thyroid function and prevent the development of thyroid diseases.Future clinical trials should investigate the therapeutic potential of vitamin D in thyroid disorders while properly accounting for the influence of season, sex and genetic background.
Although mechanisms by which 25(OH)D concentration might affect TSH levels are not fully understood, the study in rat pituitary cells showed that administration of calcitriol affects thyrotropin-releasing hormone (TRH)-induced TSH release [41].Additionally, studies in rat thyroid cells [42] and FRTL-5 cells [43] showed that calcitriol suppresses TSH-stimulated adenylyl cyclase activity and iodide uptake, respectively.In vivo studies showed that vitamin D affects Dio2 expression.This enzyme, which is necessary for the conversion of T4 into T3, showed an increase in the expression upon vitamin D3 administration in diabetic rats.This was followed by an increase in fT3 levels and a decrease in fT4 levels [44].However, vitamin D receptor (VDR) knockout mice had only moderately reduced TSH levels, with thyroid physiology not being significantly affected [45].
To date, only two studies using MR methodology investigated the association between serum 25(OH)D concentration and thyroid function.The study conducted by Chen et al. with 10,636 participants from China observed a negative causal relationship between serum 25 (OH)D concentration and levels of TPOAb [20].This causal relationship was not bi-directional since genetically predicted TPOAb levels were not causally affecting serum 25(OH)D concentration [20].Another study, by Ye et al. in 326,409 Europeans from the UK Biobank, analyzed the causal relationship between serum 25(OH)D concentration and 106 diseases/ traits.Among the disorders examined, they found no link between the risk of autoimmune thyroid diseases (hyperthyroidism and hypothyroidism) or thyroid cancer and the genetically predicted serum 25(OH)D concentration using MR technology [19].Another recent study, published as a preprint [21] has investigated the causal association between serum 25(OH)D concentration and reference-range TSH, fT4 and the risk of hypothyroidism and Hashimoto's thyroiditis.However, the authors used an earlier GWAS on TSH and fT4 [46] which had 54,288 participants in the TSH GWAS and 49,269 in the fT4 GWAS.Here we utilize the latest published GWAS on these parameters, with 271,040 participants in TSH and 119,120 participants in the fT4 GWAS, which increases the power of the MR study.In addition, the authors [21] used GWAS for hypothyroidism and Hashimoto's thyroiditis performed in the UK Biobank and FinnGen.This can lead to biased causal estimates because of sample overlap.Ideally, the exposure and outcome datasets in a two-sample MR setting should not be overlapping.
Here we utilize the FinnGen GWAS for autoimmune hypothyroidism and hyperthyroidism to avoid sample overlap between the exposure and outcomes.
Certain vitamin D receptor (VDR) gene polymorphisms, namely TaqI (rs731236) and BsmI (rs1544410) have been found to be associated with the risk of autoimmune thyroid diseases [47].Inclusion of these polymorphisms in the MR analysis would directly violate the Exclusion assumption and would bias the results of MR, as we would have horizontal pleiotropy.In our analysis, these polymorphisms were not present in the set of instruments, therefore they do not represent a potential source of bias.
To the best of our knowledge, our study is the first MR study investigating the association between serum 25(OH)D concentration and fT3, TT3, low TSH and high TSH.
One limitation of our study is that, because of the scarcity of large GWA studies in non-European ancestries, our estimation of causal effects was limited to using only GWAS performed in individuals of European ancestry.This constraint limited our ability to estimate causal effects across diverse ethnic backgrounds.Future research should aim to examine these causal relationships in various ethnic groups.Furthermore, our analysis yielded suggestive causal results only within the analyses utilizing the Revez et al. genome-wide significant SNPs.However, our aim was to provide a comprehensive MR analysis taking into account all the available resources, therefore, we believe that the non-significant findings within the Manousaki et al. analyses will provide valuable insights to researchers investigating this topic.
Furthermore, an investigation of the sex-specific causative effects of vitamin D on thyroid function was outside the scope of our study due to the absence of sex-stratified summary statistics.
In summary, our thorough MR analysis, employing both frequentist and Bayesian methods, provided indications of a potentially beneficial causal effect of serum vitamin D on thyroid disorders.Genetically predicted serum vitamin D levels were found to have a suggestive causal effect on reducing the risk of high TSH and autoimmune hypothyroidism.

Fig 1 .
Fig 1. Scatter plots depicting the relationship of A) SNP effects on serum 25(OH)D levels against the SNP effects on high TSH; B) SNP effects on serum 25(OH)D levels against the SNP effects on autoimmune hypothyroidism.MR results were derived from the multiplicative random effects inverse variance weighted (light blue), weighted median (light green), weighted mode (dark green) and MR-Egger (dark blue) method.Each line's slope represents the estimated causal effect.TSH, thyroid-stimulating hormone.https://doi.org/10.1371/journal.pone.0304253.g001
Table in S2 Appendix).The Steiger directionality test validated the causal direction of our analysis (S8 Table in S2 Appendix).Additionally, there was no indication of heterogeneity among individual variant estimates (Cochran's Q = 118.60,p = 0.29) nor evidence of uncorrelated horizontal pleiotropy (MR-PRESSO Global test p = 0.307) (S5 Table in S2 Appendix) Table in S2 Appendix).The Steiger directionality test affirmed the causal direction of our analysis (S8 Table in S2 Appendix).